This markdown provides comparison between three meQTL analysis. All meQTL analysis were done using the R package MatrixEQTL.
For all analysis, the same DNAm data were used which underwent following QC steps:
DNAm QC
The genotype data were different for each analysis:
The SNP data were used where the uncertain genotype was replaced with the most likely genotype
The uncertain genotypes were left as they are
SNP probabilities, i.e. dosages, were used
For all three analysis, the imputed QC SNP data (no LD pruning, the MHC regions are included) were used. The data includes 3’957’338 SNPs.
All meQTLs passed the significance threshold of FDR = 0.05.
x %>%
ggplot(aes(y = cnt, x = Var2, fill = as.factor(as.numeric(Var2)))) + geom_col() + geom_text(aes(label = comma(cnt,
accuracy = 1L), y = cnt), stat = "identity", vjust = -0.2, size = 2.5) + scale_y_continuous(labels = scientific) +
# scale_fill_viridis(discrete = TRUE, alpha = 0.9, option = 'C') +
labs(title = " ", y = "Count", x = " ") + facet_wrap(~Var1, ncol = 3) + theme(legend.position = "none",
legend.title = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(),
plot.title = element_text(size = 10), axis.title = element_text(size = 8), axis.title.x = element_blank()) +
scale_fill_manual(values = cbPalette[c(1, 2)])
fdr.thr <- 0.05
plot.title <- paste0("Significant at FDR <= ", fdr.thr, " meQTLs.", "\nResult of MatrixEQTL")
GetScatterPlot2(meqtl.all.prob.df, plot.title = plot.title, cbPalette = cbPalette)
meqtls.meqtls <- list(delta = meqtl.all.prob.df[treatment == "delta", meQTL_ID], veh = meqtl.all.prob.df[treatment ==
"baseline", meQTL_ID])
intersect.meqtls <- intersect(meqtls.meqtls$delta, meqtls.meqtls$veh)
perc.olap.delta.with.veh <- scales::percent(length(intersect.meqtls)/length(meqtls.meqtls$delta), accuracy = 0.1)
perc.olap.veh.with.delta <- scales::percent(length(intersect.meqtls)/length(meqtls.meqtls$veh), accuracy = 0.1)
plot.title <- "Number of intersections between delta and baseline meQTLs"
ggVennDiagram::ggVennDiagram(meqtls.meqtls, category.names = c(paste0("GR-response, ", perc.olap.delta.with.veh),
paste0("Baseline, ", perc.olap.veh.with.delta)), label_alpha = 0.7, edge_size = 0, set_geom = "text",
set_color = "black", label = "count") + theme(legend.position = "none", legend.title = element_blank(),
panel.background = element_blank(), plot.title = element_text(size = 10, face = "italic")) + labs(x = "",
y = "", title = plot.title) + scale_color_manual(values = cbPalette) + scale_fill_gradient(low = alpha(cbPalette[2],
0.5), high = alpha(cbPalette[1], 0.5))
meqtls.snps <- list(delta = meqtl.all.prob.df[treatment == "delta", SNP], veh = meqtl.all.prob.df[treatment ==
"baseline", SNP])
intersect <- intersect(meqtls.snps$delta, meqtls.snps$veh)
perc.olap.delta.with.veh <- scales::percent(length(intersect)/length(meqtls.snps$delta), accuracy = 0.1)
perc.olap.veh.with.delta <- scales::percent(length(intersect)/length(meqtls.snps$veh), accuracy = 0.1)
plot.title <- "Number of intersections between delta and baseline meQTL SNPs"
ggVennDiagram::ggVennDiagram(meqtls.snps, category.names = c(paste0("GR-response, ", perc.olap.delta.with.veh),
paste0("Baseline, ", perc.olap.veh.with.delta)), label_alpha = 0.7, edge_size = 0, set_geom = "text",
set_color = "black", label = "count") + theme(legend.position = "none", legend.title = element_blank(),
panel.background = element_blank(), plot.title = element_text(size = 10, face = "italic")) + labs(x = "",
y = "", title = plot.title) + scale_color_manual(values = cbPalette) + scale_fill_gradient(low = alpha(cbPalette[2],
0.5), high = alpha(cbPalette[1], 0.5))
meqtls.cpg <- list(delta = meqtl.all.prob.df[treatment == "delta", CpG_ID], veh = meqtl.all.prob.df[treatment ==
"baseline", CpG_ID])
intersect <- intersect(meqtls.cpg$delta, meqtls.cpg$veh)
perc.olap.delta.with.veh <- scales::percent(length(intersect)/length(meqtls.cpg$delta), accuracy = 0.1)
perc.olap.veh.with.delta <- scales::percent(length(intersect)/length(meqtls.cpg$veh), accuracy = 0.1)
plot.title <- "Number of intersections between delta and baseline meQTL CpGs"
ggVennDiagram::ggVennDiagram(meqtls.cpg, category.names = c(paste0("GR-response, ", perc.olap.delta.with.veh),
paste0("Baseline, ", perc.olap.veh.with.delta)), label_alpha = 0.7, edge_size = 0, set_geom = "text",
set_color = "black", label = "count") + theme(legend.position = "none", legend.title = element_blank(),
panel.background = element_blank(), plot.title = element_text(size = 10, face = "italic")) + labs(x = "",
y = "", title = plot.title) + scale_color_manual(values = cbPalette) + scale_fill_gradient(low = alpha(cbPalette[2],
0.5), high = alpha(cbPalette[1], 0.5))
x %>%
ggplot(aes(y = cnt, x = Var2, fill = as.factor(as.numeric(Var2)))) + geom_col() + geom_text(aes(label = comma(cnt,
accuracy = 1L), y = cnt), stat = "identity", vjust = -0.2, size = 2.5) + scale_y_continuous(labels = scientific) +
# scale_fill_viridis(discrete = TRUE, alpha = 0.9, option = 'C') +
labs(title = " ", y = "Count", x = " ") + facet_wrap(~Var1, ncol = 3) + theme(legend.position = "none",
legend.title = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(),
plot.title = element_text(size = 10), axis.title = element_text(size = 8), axis.title.x = element_blank()) +
scale_fill_manual(values = cbPalette[c(1, 2)])
x %>%
ggplot(aes(y = cnt, x = Var2, fill = as.factor(as.numeric(Var2)))) + geom_col() + geom_text(aes(label = comma(cnt,
accuracy = 1L), y = cnt), stat = "identity", vjust = -0.2, size = 2.5) + scale_y_continuous(labels = scientific) +
# scale_fill_viridis(discrete = TRUE, alpha = 0.9, option = 'C') +
labs(title = " ", y = "Count", x = " ") + facet_wrap(~Var1, ncol = 3) + theme(legend.position = "none",
legend.title = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(),
plot.title = element_text(size = 10), axis.title = element_text(size = 8), axis.title.x = element_blank()) +
scale_fill_manual(values = cbPalette[c(1, 2)])
meqtls.delta <- list(prob = meqtl.all.prob.df[treatment == "delta", meQTL_ID], miss = meqtl.all.na.df[treatment ==
"delta", meQTL_ID])
meqtls.veh <- list(prob = meqtl.all.prob.df[treatment == "baseline", meQTL_ID], miss = meqtl.all.na.df[treatment ==
"baseline", meQTL_ID])
euler.tbl <- euler(meqtls.delta)
intersect <- euler.tbl$original.values["prob&miss"]
perc.olap.prob.with.miss <- scales::percent(intersect/length(meqtls.delta$prob), accuracy = 0.1)
perc.olap.miss.with.prob <- scales::percent(intersect/length(meqtls.delta$miss), accuracy = 0.1)
plot(euler.tbl, quantities = list(type = c("counts", "percent")), fills = c("lightgrey", "red", "brown"),
labels = c("probabilities", "missigness"))
euler.tbl <- euler(meqtls.veh)
intersect <- euler.tbl$original.values["prob&miss"]
perc.olap.prob.with.miss <- scales::percent(intersect/length(meqtls.veh$prob), accuracy = 0.1)
perc.olap.miss.with.prob <- scales::percent(intersect/length(meqtls.veh$miss), accuracy = 0.1)
plot(euler.tbl, quantities = list(type = c("counts", "percent")), fills = c("lightgrey", "brown"), labels = c("probabilities",
"missigness"))
meqtls.delta <- list(prob = meqtl.all.prob.df[treatment == "delta", meQTL_ID], mode = meqtl.all.mode.df[treatment ==
"delta", meQTL_ID])
meqtls.veh <- list(prob = meqtl.all.prob.df[treatment == "baseline", meQTL_ID], mode = meqtl.all.mode.df[treatment ==
"baseline", meQTL_ID])
euler.tbl <- euler(meqtls.delta)
plot(euler.tbl, quantities = list(type = c("counts", "percent")), fills = c("lightgrey", "brown"), labels = c("probabilities",
" most likely genotype"))
meqtls.delta <- list(prob = meqtl.all.prob.df[treatment == "delta", meQTL_ID], miss = meqtl.all.na.df[treatment ==
"delta", meQTL_ID], mode = meqtl.all.mode.df[treatment == "delta", meQTL_ID])
unique.mode.meqtl <- setdiff(meqtls.delta$mode, union(meqtls.delta$prob, meqtls.delta$miss))
unique.mode.meqtl.df <- meqtl.all.mode.df[treatment == "delta"][meQTL_ID %in% unique.mode.meqtl]
selected.meqtl <- unique.mode.meqtl.df[1]
kable(selected.meqtl %>%
select(SNP, CpG_ID, FC = beta, `t-stat`, `p-value`, FDR = fdr, Treatment = treatment) %>%
mutate_if(is.numeric, funs(as.character(signif(., 3)))))
| SNP | CpG_ID | FC | t-stat | p-value | FDR | Treatment |
|---|---|---|---|---|---|---|
| rs9961090 | cg20071505 | 0.04 | 9.99 | 3.19e-19 | 2.98e-13 | delta |
ProcessGetBoxPlot(methyl.beta.veh.df, methyl.beta.dex.df, snp.df, selected.meqtl, fdr.thr = 0.05)
unique.meqtl <- setdiff(meqtls.delta$prob, union(meqtls.delta$mode, meqtls.delta$miss))
unique.meqtl.df <- meqtl.all.prob.df[treatment == "delta"][meQTL_ID %in% unique.meqtl]
selected.meqtl <- unique.meqtl.df[1]
kable(selected.meqtl %>%
select(SNP, CpG_ID, FC = beta, `t-stat`, `p-value`, FDR = fdr, Treatment = treatment) %>%
mutate_if(is.numeric, funs(as.character(signif(., 3)))))
| SNP | CpG_ID | FC | t-stat | p-value | FDR | Treatment |
|---|---|---|---|---|---|---|
| rs11655681 | cg17105886 | -0.0407 | -5.75 | 3.35e-08 | 0.00231 | delta |
ProcessGetBoxPlot(methyl.beta.veh.df, methyl.beta.dex.df, snp.df, selected.meqtl, fdr.thr = 0.05)
# meqtls.delta <- list(prob = meqtl.all.prob.df[treatment == 'delta', meQTL_ID], miss =
# meqtl.all.na.df[treatment == 'delta', meQTL_ID], mode = meqtl.all.mode.df[treatment == 'delta',
# meQTL_ID])
meqtl.lst <- intersect(intersect(meqtls.delta$mode, meqtls.delta$prob), meqtls.delta$miss)
meqtl.df <- meqtl.all.mode.df[treatment == "delta"][meQTL_ID %in% meqtl.lst]
selected.meqtl <- meqtl.df[1]
kable(selected.meqtl %>%
select(SNP, CpG_ID, FC = beta, `t-stat`, `p-value`, FDR = fdr, Treatment = treatment) %>%
mutate_if(is.numeric, funs(as.character(signif(., 3)))))
| SNP | CpG_ID | FC | t-stat | p-value | FDR | Treatment |
|---|---|---|---|---|---|---|
| rs9406094 | cg12777182 | 0.0694 | 16.6 | 3.06e-39 | 3.01e-30 | delta |
selected.meqtl <- meqtl.df[1]
ProcessGetBoxPlot(methyl.beta.veh.df, methyl.beta.dex.df, snp.df, selected.meqtl, fdr.thr = 0.05)
meqtl.df <- meqtl.all.prob.df[treatment == "delta"][meQTL_ID %in% meqtl.lst]
selected.meqtl <- meqtl.df[1]
kable(selected.meqtl %>%
select(SNP, CpG_ID, FC = beta, `t-stat`, `p-value`, FDR = fdr, Treatment = treatment) %>%
mutate_if(is.numeric, funs(as.character(signif(., 3)))))
| SNP | CpG_ID | FC | t-stat | p-value | FDR | Treatment |
|---|---|---|---|---|---|---|
| rs12310416 | cg07195891 | -0.0718 | -16.7 | 2.35e-39 | 1.84e-30 | delta |
selected.meqtl <- meqtl.df[1]
ProcessGetBoxPlot(methyl.beta.veh.df, methyl.beta.dex.df, snp.df, selected.meqtl, fdr.thr = 0.05)
meqtl.df <- meqtl.all.na.df[treatment == "delta"][meQTL_ID %in% meqtl.lst]
selected.meqtl <- meqtl.df[1]
kable(selected.meqtl %>%
select(SNP, CpG_ID, FC = beta, `t-stat`, `p-value`, FDR = fdr, Treatment = treatment) %>%
mutate_if(is.numeric, funs(as.character(signif(., 3)))))
| SNP | CpG_ID | FC | t-stat | p-value | FDR | Treatment |
|---|---|---|---|---|---|---|
| rs12310416 | cg07195891 | 0.0718 | 16.7 | 2.35e-39 | 3.21e-30 | delta |
selected.meqtl <- meqtl.df[1]
ProcessGetBoxPlot(methyl.beta.veh.df, methyl.beta.dex.df, snp.df, selected.meqtl, fdr.thr = 0.05)